Quantum quenches and driven dynamics in a single-molecule device 
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The nonequihbrium dynamics of molecular devices is studied in the hamework of a generic model 
for single-molecule transistors: a resonant level coupled by displacement to a single vibrational 
mode. In the limit of a broad level and in the vicinity of the resonance, the model can be controUably 
reduced to a form quadratic in bosonic operators, which in turn is exactly solvable. The response of 
the system to a broad class of sudden quenches and ac drives is thus computed in a nonperturbative 
manner, providing an asymptotically exact solution in the limit of weak electron-phonon coupling. 
From the analytic solution we are able to (1) explicitly show that the system thermalizes following 
a local quantum quench, (2) analyze in detail the time scales involved, (3) show that the relaxation 
time in response to a quantum quench depends on the observable in question, and (4) reveal how 
the amplitude of long-time oscillations evolves as the frequency of an ac drive is tuned across the 
resonance frequency. Explicit analytical expressions are given for all physical quantities and all 
nonequihbrium scenarios under study. 

PACS numbers: 73.63.b, 71.38.k, 85.65.-(-h 
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I. INTRODUCTION 

The description of strong electronic correlations far 
from thermal equilibrium constitutes one of the ma- 
jor open questions of modern condensed matter physics. 
Even under the most favorable conditions of nonequilib- 
rium steady state, many of the concepts and techniques 
that have proven so successful in equilibrium are sim- 
ply inadequate. Recent advancements in a broad range 
of systems, from time-resolved spectroscopiesiS to cold 
atoms^il and driven nanostructures?^ have opened new 
and exciting possibilities for studying the nonequilibrium 
dynamics in response to quantum quenches and forcing 
fields. Depending on the physical context one is inter- 
ested in questions of both basic and practical nature, such 
as what are the underlying time scales governing the dy- 
namics, how long is coherence maintained, and whether 
and how does the system equilibrate at long times. Some 
questions, e.g., the issue of equilibration, often require 
nonperturbative treatments even if the system is tuned 
to weak coupling. 

Recent years have witnessed the development of an 
array of powerful numerical techniques aimed at track- 
ing the real-time dynamics of interacting low-dimensional 
systems. In the more specific context of quantum impu- 
rity systems these methodologies include time-dependent 
variants of the density-matrix renormalization group, ^'^ 
the time-dependent numerical renormalization group^iiS 
different continuous-time Monte Carlo approaches r^"— 
and sparse polynomial space representations^^ Despite 
notable successes, part of these methods are subject to 
finite-size effects and discretization errors, while others 
are confined to rather short time scales. Analytical ef- 
forts in this realm have focused mainly on suitable adap- 
tations of perturbative renormalization-group^^*i^ and 
flow-equationi^ ideas, which in turn neglect higher or- 
der terms. Exact analytical solutions, when available, 



are thus invaluable both for setting a benchmark and 
for gaining unbiased understanding of the underlying 
physics. Unfortunately such exact solutions are restricted 
at present to very special models whose coupling con- 
stants must be carefully tuned.— 

In this paper we present an asymptotically exact solu- 
tion for the nonequilibrium dynamics of a single-molecule 
transistor in response to various quantum quenches and 
ac drives. Single-molecule devices have attracted consid- 
erable interest lately due to the technological promise of 
molecular electronics.^^ From a basic-science perspective 
they offer an outstanding platform to study the electron- 
phonon coupling at the nano-scale. In a typical molecular 
bridge, molecular orbitals are coupled simultaneously to 
the lead electrons and to the vibrational modes of the 
molecule, with the former degrees of freedom reduced 
to a single effective band in the absence of a bias volt- 
age. A minimal model for an unbiased molecular bridge 
therefore consists of a single resonant level coupled by 
displacement to a single vibrational mode, as described 
by the Hamiltonian of Eqs. (H]) and ^ below. 

The spinless Hamiltonian of Eqs. ([T]) and ^ has been 
extensively used in recent years to model single-molecule 
transistors, however despite its apparent simplicity it 
lacks a complete solution. Conventionally the model is 
treated either using perturbation theory in the electron- 
phonon coupling when the coupling is sufficiently weak, 
or using the Lang-Firsov transformation'^"^ and the po- 
laronic approximation in the limit where tunneling is 
sufficiently small. A particularly elegant nonperturba- 
tive solution of the model was recently devised by Dora 
and Halbritter,^"* who noticed that the original electronic 
Hamiltonian of Eqs. ^ and ([2]) can be mapped onto 
an exactly solvable bosonic form in the limit where the 
electronic level is broad. Building on prior result o^^'^^ 
for the related single-impurity Holstein model, these au- 
thors proceeded to compute the temperature-dependent 
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conductance of the device under strict resonance condi- 
tions. Since mapping onto the exactly solvable model is 
controlled by the smallness of the electron-phonon cou- 
pling g as compared to the level width F, these results are 
expected to be asymptotically exact in the weak-coupling 
limit. 

In this paper we take the solution one step further 
by extending it to the nonequilibrium dynamics in re- 
sponse to a broad class of quantum quenches and drives. 
We explicitly show that the system thermalizes following 
a local quantum quench and analyze in detail the time 
scales involved. In particular, we find that the relaxation 
time depends on the observable in question, growing by a 
factor of two in going from the phonon occupancy to the 
phonon displacement and the electronic occupancy of the 
level. This is quite surprising since unlike the Anderson 
impurity model, where spin and charge generally relax 
on different time scales^ the phonon occupancy and dis- 
placement pertain to the same degrees of freedom. A 
related doubling of frequency occurs in the long-time re- 
sponse of the phonon occupancy to an ac drive. These 
results, as well as others, are obtained in a fully analytic 
manner, which is perhaps the most appealing aspect of 
our solution. 

Before proceeding to actual calculations, two technical 
comments are in order. First, some of the scenarios under 
consideration in this paper pertain to a level off resonance 
with the Fermi energy, which necessitates the incorpora- 
tion of the level energy into the bosonic Hamiltonian. 
A nonzero energy level breaks particle-hole symmetry, 
an aspect that is missing in the treatment of Dora and 
Halbritter. Below we correct their mapping to properly 
account for this important point. Second, some of the 
initial states to be considered will be nonthermal states 
that cannot be treated using the Keldysh technique. We 
circumvent this complication by explicitly constructing 
the single-particle eigenmodes of the bosonic Hamilto- 
nian and using them to propagate the system in time. 

The reminder of the paper is organized as follows. In 
Sec. HI] we introduce the model and its mapping onto a 
form quadratic in bosonic operators. The bosonic Hamil- 
tonian is solved in turn in Sec. lIIII bv explicitly construct- 
ing its single-particle eigenmodes using the scattering- 
state formalism. Technical details of the solution are 
relegated to the Appendix. The next three sections 
are devoted to three different quench scenarios: one. 
Sec. IIVI where the electron-phonon interaction is sud- 
denly switched on, another. Sec. |Vl where the phonon 
frequency is abruptly shifted from its initial value, and 
lastly. Sec I VII the scenario where a sudden change is ap- 
plied to the electronic level. The case of driven dynamics 
is addressed in Sec. IVIIl first in its general form before 
turning to an explicit discussion of ac drives. Finally, we 
present our conclusions in Sec. I Villi 



II. THE MODEL AND ITS MAPPING 

The Hamiltonian we consider is one of the common 
models used to describe a single Coulomb-blockade reso- 
nance in molecular devices. It consists of a single spinless 
electronic level with energy e^, which is coupled by dis- 
placement to a local vibrational mode with frequency 
wq. The level is further coupled to a band of spinless 
electrons via the hopping matrix element t, as described 
by the Hamiltonian^! 

n = Ho + edfid + oj^b^b + g {b^ + b) (^a - ^ , (1) 
with hd = d^d and 

no^J2 ^kcick + t {eld + d^Cfc) . (2) 

k k 

Here the combination Q = (6^ + fe)/\/2 can be thought of 
as a dimensionless position operator for the local phonon. 

The Hamiltonian defined by Eqs. ([T]) and ([U has a 
long history that dates back to the 1970s, when it was 
proposed as a model for the electron-phonon coupling 
in mixed-valence compoundsj^ In the modern context 
of nanostructures it is expected to properly describe the 
physics of single-molecule devices away from Coulomb- 
blockade valleys where a single unpaired spin resides on 
the molecule. Typically the Hamiltonian is treated either 
in the weak-coupling limit using perturbation theory in 
g, or using the Lang-Firsov transformation^'^ and the po- 
laronic approximation in the limit where t is small. We 
shall take a different route and present a nonperturbative 
solution to this model which is asymptotically exact in 
the limit where F 3> maxjg, |erf|, g^/wo}. Our approach 
is based on the fact that the Hamiltonian of Eqs. ([T]) and 
can be controUably reduced in this limit to a form 
quadratic in bosonic operators which is exactly solvable. 
As discussed in the introduction, this method was first 
employed in equilibrium by Dora and Halbritter Here 
we exploit this property of the model to calculate the 
real-time dynamics following different quantum quenches 
and also in response to ac drives. Accordingly, our pre- 
sentation begins with the conversion of the Hamiltonian 
to a form that is quadratic in bosonic operators, whose 
solution is detailed in turn in Sec. IIIII 

Technically, the construction of the bosonic Hamilto- 
nian proceeds in two steps (Ref . (H) : (i) the conversion 
to a continuum-limit Hamiltonian and (ii) its subsequent 
bosonization. Special care is paid to the parametric form 
of the coupling constants that enter the bosonic Hamil- 
tonian and to the role of the energy level td which breaks 
particle-hole symmetry. The latter energy scale is of par- 
ticular interest as it can be controlled experimentally us- 
ing suitable gate voltages. In this respect our derivation 
exceeds that of Dora and Halbritter. 
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A. Conversion to a continuum-limit Hamiltonian 

Our first goal is to map the Hamiltonian of Eq. ([ij onto 
a continuum-limit form, where the resonant-level opera- 
tor is replaced with a suitable field operator. To this 
end, we first diagonalize the Hamiltonian term Ho using 
scattering theory to construct its single-particle eigen- 
modes. These are conveniently expressed using the Green 
function of the level 



G(z) 



and its associated phases 

arg{G(efe - ir])} , 



(3) 



(4) 



where the limit 77 — > 0+ is implied. Specifically, introduc- 
ing the properly normalized fcrmionic operators 



d^ + y 



(5) 



the Hamiltonian term H.^ can be shown to take the diag- 
onal form 



(6) 



Our manipulations thus far were exact, independent 
of details of the band dispersion tk- To make further 
progress we consider hereafter the wide-band limit, where 
the spectral function of Eq. (|lip acquires the Lorentzian 
form Trpd{e) = r/(e^ + F^) with the hybridization width 
r = TTp{0)P [p{0) is the conduction electrons density of 
states at the Fermi energy] . Physically, F serves as a new 
high-energy cutoff for the integration in Eq. ([TU]). Since 
d^ is the only electronic degree of freedom that enters 
the remaining Hamiltonian terms in Eq. ([1]), F acts as 
a new effective bandwidth for the electron-phonon cou- 
pling. We shall next exploit this observation to further 
manipulate the Hamiltonian of the system. 

The Lorentzian cutoff in Eq. ([TU]) is somewhat inconve- 
nient to deal with. However, its precise form should not 
play any role in the desired limit F 3> max{(7, \ ed\,g'^/uJo}, 
allowing one to adopt a more convenient cutoff scheme. 
Indeed, it is useful to replace Pd{^) in Eq. (ITUl) with a rect- 
angular box profilc^^ that has the same height at e = 
and shares the same characteristic width: 



with 



i{e) ^ j^9{Dd 



ttF 

Dd = ^. 



(12) 



(13) 



Substituting P(j(e) with the box profile of Eq. (fT2)) . the 
full Hamiltonian of Eq. ([1]) becomes 



while d^ acquires the mode expansion 

d^ ^tJ2\Giek+^v)\4■ 



(7) 



Further converting to the continuous energy-shell opera- 
tors 



^P(^) k 



(8) 



where p(e) is the conduction-electron density of states, 
Eqs. ([6]) and (O become 



etpltp^de 



and 



d^ 



D 



^/ Pd{e)tpld€. 



Here D is the conduction-electron bandwidth and 



Pd{e) = Im{G(efc + irj)} 

TT 



(9) 



(10) 



(11) 



is the spectral function associated with the Green func- 
tion of Eq. ©. 
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(14) 



de de' -.i^li,,,:^ 



where :'4'l'4'e'-= V'lV'e' ~ ^{^ ~ £')^(~£) stands for normal 
ordering with respect to the filled Fermi sea. Note that 
all electronic modes with |e| > Dd are decoupled from the 
phonon in Eq. (|14p and can therefore be omitted. This 
amounts to setting D — >■ Dd in the integration boundaries 
of the free kinetic-energy term. 

The conversion to a continuum-limit Hamiltonian is 
completed by defining the right-moving field 



1 



V^MDd J -Da 



where is the Fermi velocity and 



■KVp 

a7 



2vf 



(15) 



(16) 



is a new short-distance cutoff corresponding to a lattice 
spacing. The new cutoff is connected to the momen- 
tum cutoff kc = vp/Dd through the standard relation 
kc — -K /a. The field operators so defined obey canonical 
anticommutation relations {^/'(a;), ^^(y)} = S{x — y), sub- 
ject to the regularization 5{0) = 1/a. Recalling that the 



4 



local fermion d'' has been mapped in this process onto 
Y^V^(O), this regularization guarantees that {d,d^} = 1 
is preserved. Written in terms of the new field operators, 
the Hamiltonian of the system takes the continuum-limit 
form 

/CO 
if)^ {x)dx'4^{x)dx + LOoh^b 
-co 

+ [e~rf + A(6t + 6)]:^t(o)^(0): (17) 

with 

A = 5" = 2^5, (18) 

id = e^a = 2— erf- (19) 

Hence, the resonance width F, the electron-phonon cou- 
pling and the energy level have been reduced to 
two parameters only, which have the dimension of en- 
ergy times length. It should be stressed that the original 
conduction-electron bandwidth D has been replaced in 
Eq. P7|) with Dd ^ F, which serves as the new high- 
energy cutoff for the continuum-limit Hamiltonian. 

The Hamiltonian of Eq. (IT7| , first derived in this con- 
text by Dora and HalbritterjSl is by no means new. It de- 
scribes the coupling of a localized phonon mode to a con- 
duction band, and as such has been applied in different 
variants to a broad class of physical systems. For exam- 
ple, Gadzuk considered it as a general impurity model'^'^ 
before applying it to the vibrational line shape of di- 
atomic adsorbates on metallic clusters.— Yu and Ander- 
son^-» proposed a closely related two-band Hamiltonian 
as a model for the anomalous properties of A15 mate- 
rials, while Dora and Gulacsi^ used this Hamiltonian 
to study the inelastic scattering from local vibrational 
modes. Although the model in its general form lacks a 
full solution, it can be conveniently handled in the pa- 
rameter regime of interest to us using the methodology 
of Abelian bosonization. 



number operator, : O : stands for normal ordering with 
respect to the filled Fermi sea, and ^ is a phase operator 
conjugate to N . The coefficients have the explicit form 

e, = /^e-/-, (22) 

which includes a suitable ultraviolet momentum cutoff 
kc = 7r/a. 

The rules of bosonization enable one to represent 
fermionic operators in terms of bosonic ones with an 
important caveat: the bosonized form of the interac- 
tion term is generally not known away from weak cou- 
pling. This uncertainty is removed in the limit of interest 
F 3> max{(7, le^l, g^/wo}, when the standard substitution 
: iIj\x)^{x) :— {~l/2T:)dx4>{x) applies. Restricting at- 
tention to this regime, the bosonized Hamiltonian is thus 
recast as^'* 

% = ekajak + ujpb^b 

k>0 

+ [A(&t+6) + f:rf]^e,K + 4)- (23) 

g>0 

Another important identity pertains to the occupancy 
of the localized electronic level fit; — d^d. Since 
has been mapped in the continuum limit onto •\/aV'^(0)j 
then fid — 1/2 corresponds to a: ■(/;^(0)V'(0) :, where we 
have made use of the fact that the expectation value of 
"0^ (0)^0(0) with respect to the unperturbed Fermi sea is 
l/(2a) [see Eq. with x = 0]. Accordingly, 7x^-1/2 
has the bosonized representation^ 

f^d~]^^aY^^k{al^ak). (24) 

This identity will play a key role in our calculations be- 
low. 

III. EXACT DIAGONALIZATION 



B. Abelian bosonization 

Our next step is to bosonize the continuum-limit 
Hamiltonian defined by Eq. (|T7|) . Using the standard pre- 
scriptions of Abelian bosonizatiouj^ the fermionic field 
operator ip{x) is written as 

V(x) = ^e-*(-), (20) 
where the bosonic field 0(a;) has the mode expansion 
</)(x) = 2^z y ^ (a„e*«^ - aje"'"'") - — : TV: + 0. 

(21) 

Here, and are canonical bosonic creation and anni- 
hilation operators corresponding to the Fourier compo- 
nents of the electronic density, N is the total fermionic 



The Hamiltonian of Eq. ([25)) is quadratic in bosonic 
operators and as a result is exactly solvable. In the fol- 
lowing section we construct its single-particle eigenmodes 
using the scattering-state formalism. Although of similar 
technical complexity, it is advantageous to first address 
the case where = 0, and then extend the discussion to 
nonzero Cd- This will prove beneficial as we shall be in- 
terested, among other things, in cases where the level en- 
ergy is shifted abruptly from = to nonzero e^. As we 
shall see, such a scenario requires the conversion between 
the eigenmodes of the Hamiltonian with and without ed- 
In contrast to the Keldysh technique, the expansion in 
terms of the eigenmodes of the bosonic Hamiltonian will 
enable us to address cases of practical interest where the 
system is initially prepared in a nonthermal state. For 
example, if the phonon initially occupies an excited state. 
Our approach is therefore more general than the Keldysh 
technique. 
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A. Scattering states for ed — 

When a free bosonic mode impinges upon the local 
phonon , it is scattered into a linear combination of the 
free bosonic modes and the localized phonon. This pro- 
cess can be described by the scattering states a^, which 
are eigenmodes of the bosonic Hamiltonian obeying suit- 
able boundary conditions of an incoming free particle. 
The scattering-state operators can be found by solving 
the Lippmann-Schwinger equation, which takes the op- 
erator form 



-ekal 



a 



(25) 



The role of — >■ 0+ in Eq. (|25p is to guarantee appropri- 
ate boundary conditions. It does not enter any physical 
quantities. 

A detailed solution of Eq. ([25]) is presented in Ap- 
pendix 13 using the methodology developed in Ref. |3^. 
Here we quote only the end result. The scattering-state 
operators are given by 




(efe - uJo)b + (ffe + wo)&^ 



— lOq — 2a;oS(z) 



(26) 



(27) 



is related to the phononic Green function of Eq. (|A20[) 
and 



E(z) = A- 



fc>0 



1 



1 



(28) 



is the corresponding self-energy. Both S(z) and g{z) are 
analytic in the upper and lower halves of the complex 
plane, have a branch cut along the real axis, and are even 
functions of z [i.e., g{z) = g{—z) and likewise for S(z)]. 
In addition g{z*) = g*{z) and E(z*) = E*(z). These 
analytical properties are useful in establishing some of 
the operator identities that will be employed in this pa- 
per. In particular, it can be explicitly shown that, in 
the limit where i — > oo, rj ^ and yet Lrj — ^ cx), the 
Hamiltonian takes the diagonal form 



(29) 



k>0 



while the scattering-state operators maintain canonical 
commutation relations: 



(30) 



In fact, the latter commutation relations apply to any 
finite ry, though only the limit ?7 — >■ 0^ is of interest to us. 

One particularly useful identity is the expansion of the 
local phonon mode in terms of the scattering-state 
operators; 



k>0 



-g{ek + j?/)(efc - ii^o)afc 



(32) 



Combined with the diagonal form of the Hamiltonian of 
Eq. (j29l) . one can immediately write down the time evo- 
lution of (t) in the Heisenberg representation, which 
reads 



k>0 



g{ek - ivi){tk+i^o)e""''^al 



-g{ek + iri){ek-uJo)e "'=*afc 



(33) 



A similar identity applies to the occupancy of the lo- 
calized electronic level, whose bosonized form has been 
detailed in Eq. (IM)) . Expanding the right-hand side of 
Eq. dUl) as 

a^£.k{.tl- ^l) g{ek - + g{ek + iv)ak , (34) 

k>0 

one has that 



k>0 



g{ek - iv)e'"''al 



+g{ek + i'n)e 



Ctk 



(35) 



The operator identities listed in Eqs. ([55)1 and ([55)1 
are central to our study as they allow one to track the 
nonequilibrium dynamics of the phonon mode and the 
level occupancy, respectively. Accordingly, they will be 
heavily used throughout the paper. 



B. Extension to nonzero td 

As stated above, the inclusion of a nonzero is quite 
straightforward and does not add to the complexity of 
computing the scattering-state operators. Since adds 
a term linear in bosonic operators to the Hamiltonian [see 
Eq. the resulting scattering-state operators differ 

by a simple fc-dependent displacement from their e^j = 
counterparts (see Appendix El for a detailed derivation) . 
Reserving the notation a\ for the scattering-state opera- 
tors when Erf = and denoting the new operators by /3]^, 
the latter are given by 



[ak.aq] = 0. 



(31) 



Pi 



al+ed^k^^—^g{ek + iri), (36) 
ek+iri 
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where and a], are specified in Eqs. pJ|) and ((26t . re- 
spectively. The shift in scattering-state operators carries 
over to physical observables as well. For example, the 
the local phonon mode is expanded as 



(37) 



where b^ is given by the same formal expression of 
Eq. ([5^ . but with al and ak replaced with and /3k, 
respectively: 



k>0 



-g{ek + ivi){ek - t^ojfik 



(38) 



As discussed below [see Eq. ([5^ with ^ — >• 0], the 
self-energy I](— ijy) takes the particularly compact form 
— (7^/(7rr), hence Eq. ([57]) can be rewritten as 



b^ = b^ + 



g 



ttE Wo 1 - 2g2/(7ra;or) ' 



(39) 



where we have expressed the constant shift in terms of 
the original model parameters that appear in Eq. ([T]). 

An analogous expansion applies to the occupancy of 
the localized level, fid, which is written as 

fid ^ fid - '^u)lg{-irii)T,{-ini). (40) 



Here 



rid 



k>0 



(41) 

is the same formal expression of Eq. ([55]) with the time 



t set to zero and with aj, and ak replaced by f3j^ and 
/3k, respectively. As with 6^, one can exploit the explicit 
expression for the self-energy I](— iry) to recast hd in the 
form 



nd = nd 



e-d 



1 



(42) 



Note that the displacement of the scattering-state op- 
erators and the associated shifts in the expansions of W 
and fid have a simple physical origin: they reflect the 
breaking of particle-hole symmetry in the original Hamil- 
tonian of Eq. ([ij inflicted by a nonzero e^. This im- 
portant aspect of is absent in the mapping of Dora 
and Halbritter^ who accounted for this energy scale by 
a simple Lorentzian reduction of the coupling constant 
A. Some of the results presented in this work would be 
missed out unless the breaking of particle-hole is properly 
treated. 

Armed with the single-particle eigenmodes of the full 
Hamiltonian and with the expansions of physical opera- 
tors in terms of these modes, we are now in position to 



compute the real-time dynamics of the system in response 
to various quantum quenches and ac drives. Specifically, 
we shall consider three quench scenarios: one where the 
electron-phonon interaction is abruptly switched on, an- 
other where the phonon frequency is suddenly shifted 
from cjQ to Wo + <5w, and finally a sudden change in the 
level energy from = to nonzero e^. In addition, we 
shall consider two ac drives — one applied to the local 
phonon and another applied to the electronic level. Of 
particular interest are the characteristic time-scales that 
govern the nonequilibrium dynamics and their depen- 
dences on the physical parameters of the system. These 
aspects will be analyzed in detail below. 



IV. SWITCHING ON THE INTERACTION 

We begin our discussion with the nonequilibrium dy- 
namics following an abrupt switching on of the electron- 
phonon interaction g. We consider the following scenario. 
At time t < the system is free of interactions (i.e., 
g = 0), and occupies a state that is a direct product of 
the electronic ground state (the filled Fermi sea) and an 
arbitrary phononic state. Typically one is interested in 
cases where the phonon has either a well-defined occu- 
pation number n or resides in a coherent state, though 
our discussion is not restricted to these particular choices. 
At time t — Q the electron-phonon interaction is abruptly 
switched on and the system evolves under the full Hamil- 
tonian T-L. This acts to entangle the phononic and elec- 
tronic degrees of freedom, which are no longer indepen- 
dent. We concentrate our discussion on zero tempera- 
ture, yet the derivation presented below can readily be 
extended to any finite temperature T of the Fermi sea. 

Formally, the time evolution of the expectation value 
of an observable O is given by the standard expression 



o{t) = {^o\u\t,Q)du{tM^o) 



(43) 



where |?/'o) is the initial state of the system and U{t, 0) is 
the time-evolution operator. One is therefore interested 
in the expectation value of O in its Heisenberg represen- 
tation d{t) = U\t, 0)dU{t, 0) with respect to the initial 
state {ipo}- In the scenario under consideration the initial 
state has a simple representation in terms of the eigen- 
modes of the initial Hamiltonian with g = 0, whereas 
the time evolution has a natural representation in terms 
of the eigenmodes of the full (i.e., final) Hamiltonian H. 
Therefore, the general strategy for calculating 0{t) pro- 
ceeds as follows. First 0{t) is represented in terms of the 
eigenmodes of the full Hamiltonian where its time evolu- 
tion can easily be implemented, next it is recast in terms 
of the eigenmodes of the initial Hamiltonian, and finally 
the expectation value with respect to \ipo) is evaluated. 
Below we implement this procedure to track the time 
evolution of the phononic occupancy n&(i) — {b^{t)b{t)) 
and displacement Q{t) = ■^{b'^ {t) + b{t)). Throughout 
the section we set Cd equal to zero, corresponding to a 
level at resonance with the Fermi energy. 
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A. Time evolution of phononic operators 

Our first goal is to express (<) in terms oiak, a\, b and 
b^ , which are the eigenmodes of the initial Hamiltonian 



with (7 = 0. The expansion of b^t) in terms of the eigen- 
modes of the final Hamiltonian is detailed in Eq. p3|) . 
Substituting the explicit expression for the scattering- 
state operators, Eq. (^5]) . into Eq. ([55]) one obtains 



b^ (i) = A ^ 6 [F{ek - ill, t)al + F{-ek - ill, t)ak\ + h {t)b + h [tp 

fc>0 

where we have introduced three auxiliary functions 

k>Q 
fe>0 

F{z, t) = g{z) [z + ^o)e'"* + ^ (I |5(e, + 

r 



A;>0 



efc -I- z 



(44) 

(45) 
(46) 

(47) 



In general, one must resort to numerical integration to 
accurately evaluate the three functions defined above at 
arbitrary time t. Results of such numerical calculations 
will be presented below for the relevant observables of 
interest. It is instructive, however, to first gain analyt- 
ical insight by analyzing the long-time behaviors of the 
auxiliary functions. In the limit L — > oo one can replace 
the sums over fc with integrals over energy, resulting in 
an exponential decay at long times of all items but the 
first term on the right-hand side of Eq. (|47)) . To see this 
important point consider I\(t\ for example. Converting 
the sum over k into integration over energy, one is left 
with the integral 



/i(i) = (poA)2 / deW^iri)\\^ -u?^) 
Jo 

xe(e"*-e-"*)e-^/^^ (48) 

where po = 1/{2t:vf) is the density of states per unit 
length. Focusing on t ^ ^/Dd, one may (i) omit the 
exponential cutoff e"^/^"* in Eq. (^51) and (ii) interchange 
e — ^ — e in the second term in the parenthesis to obtain 

/>oo 

/i(i)^(poA)M de\g{e + ^r^)\\{e'-ul)e^^K (49) 



The function g{e + it]) is analytic in the upper half of 
the complex plane, whereas the analytic continuation of 
g*{€ + irf) to the upper half plane has a set of isolated 
poles^ of the form pj = ujj + i/tj with tj > 0. Us- 
ing these poles one can formally perform the integral in 
Eq. (j49|) to arrive at 



hit) ~ 27rz(poA)2 J2 Rj{p] - ^o)Pje 



where Rj is the residue of \g{e + ir?)p at pj. Thus, for 
t ^ l/Dd, the function Ii{t) is well approximated by a 
discrete sum of exponential terms that contain both an 
oscillatory component and a part that decays in time. 
Asymptotically only those terms with the largest decay 
time Tj dominate, hence Ii(t) closely follows a simple ex- 
ponential decay with superimposed oscillations. A sim- 
ilar procedure can be applied to hit) and to the term 
involving the sum over k in the expression for i^(z,t), 
both of which arc found to be dominated by the same set 
of poles Pj provided z lies in the lower half plane (as is 
the case throughout our calculations). 

Next we address the poles pj, which are given by the 
solutions to the equation 



z^-uj^- 2woS(+Hz) = 0, 



(51) 



where S](+)(z) is the analytic continuation of E*(e -I- irj) 
to the upper half plane. For L — oo, the self-energy of 
Eq. P5|) has the explicit analytic expression 



S(z) = ipoX)'Dd [Ce«i?i iO - ^e-^E, (-^ - 2] , (52) 

where ^ equals z/Dd and Eiiz) is the Exponential 
Integral functioni^ Expanding Eiiz) as a logarithm 
plus a power series in z one obtains E*(e -I- irf) = 
ipoX)'^Dd [iTre - 2 + O (e^ Ine)] with e = e/Dd, resulting 
in 

iz) = ipoX)^Dd [ittC -2 + 0{e Ine)] • (53) 



iujjt—t/rj 



(50) 



In general, Eq. ([5T|) lacks an analytical solution. How- 
ever, in the desired limit where p^X = g/iT:T) ^ 1 and 
^ E)d one can truncate the expansion of E(+)(z) at 
linear order in e, to be left with a simple quadratic equa- 
tion in Eq. (|5ip . In this limit only two poles exist, which 
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differ in the sign preceding the frequency: p± = ztw + z/r 
with 



(54) 



The decay time r is given in this approximation by 
T = 7rr^/(cjo<?^), where we have converted back to the 
original model parameters of Eq. ([1]) in writing both the 
frequency uj and the single relaxation time t. Note that 
these expressions for w and r coincide with second-order 
perturbation theory in g when applied directly to the 
electronic Hamiltonian of Eq. thus validating the 

cutoff scheme used in bosonization. The expression for 
r can be further improved by going to the next order 
in ujo/Dd, i.e., by including one more order in f in the 
expansion of S^+'(z). This in turn yields 



^0 \g 



2 UJn 



(55) 



where we have restricted ourselves to linear order in 
uq/Dci in writing the expression in the square brackets. 

As we shall confirm by explicit numerical calculations, 
the nonequilibrium dynamics of all observables of interest 
is governed exclusively by uj and t at time scales exceed- 
ing l/Dd- Similar results for uj and l/r were reported by 
Dora and Halbritter — (corresponding in their notation to 
the real and imaginary parts of ujp± ) , yet their expression 
for u! contained the bare conduction-electron bandwidth 
D rather than the renormalized one Dd ^ T D. In- 
deed, Eqs. (IMl) and ([55)1 are free of the bandwidth D, in- 
dicating that one can safely implement the limit I? — > cx) 
for the Hamiltonian of Eq. ([T|), provided F, g, and ujq 
are all held fixed. Physically this refiects the fact that 
the local phonon couples to the conduction band by way 
of the resonant level only, hence its level width F serves 
as a new effective high-energy cutoff. By contrast, there 
is no meaningful Dd oo limit for the continuum-limit 
Hamiltonian of Eq. p7)) that keeps both ut and t finite. 

Equation ([M)) features two special values of the 
electron-phonon coupling g. One, go, above which the 
frequency uj becomes imaginary (i.e., p± become purely 
imaginary) and another, slightly larger coupling gc, above 
which the pole p- is shifted to the lower half plane. 
The former coupling strength represents the point where 
the local phonon is completely softened, whereas the 
latter value represents the point above which the en- 
ergy of the lowest bosonic eigenmode of the Hamilto- 
nian of Eq. becomes negative. Both values of g 
lie well beyond the applicability of our theory, as the 
mapping onto the continuum-limit Hamiltonian assumed 
F ^ max{g, g^ /ujq}. Interestingly, it has been argued by 
Dora^^ that the bosonized Hamiltonian of Eq. ([25)1 offers 
a faithful representation of the continuum-limit Hamilto- 
nian of Eq. (jl7|) all the way up to strong coupling, where 
the nonlinear conversion between the fermionic and the 
bosonic coupling constants is not explicitly known. In 
particular, the point where p_ is shifted to the lower 



half plane was identified by Dora with A — oo. It re- 
mains to be seen whether such strong electron-phonon 
couplings can indeed be described by a bosonized Hamil- 
tonian with just a simple linear displacement coupling, 
or whether additional nonlinear terms must be included. 



B. Phononic occupancy and displacement 

With the explicit expansions of b^{t) and b{t) at hand 
we can proceed to compute the time evolution of physical 
observables, starting with the phonon occupancy nb{t) 
and displacement Q{t). Since b^{t) is linear in the eigen- 
modes of the initial Hamiltonian, the phonon number op- 
erator hii{t) — {t)b{t) is quadratic in these operators. 
When averaged with respect to the initial state, only the 
combinations Ofcaj^, b^b, bb\ bb, and 6^6^ contribute to 
the expectation value of nf, at time resulting in 



fc>0 



+ \h{t)V nM + 2Re{h{t)I*{t){bb)t=o} ■ (56) 



rib 



Equation (156)) is the central result of this subsection. 
It provides an asymptotically exact expression for the 
time evolution of nh{t) in the weak-coupling regime. Sev- 
eral comments should be made about this result. First, 
the occupancy nb(i) depends on the initial state of the 
phonon via two parameters only: n{,{0) and {bb)t=o- Any 
two initial states that share the same values of nf,(0) and 
{bb)t=o will produce identical curves for nf,{t). Second, 
since Ii{t) and /2(i) decay to zero with time, the occu- 
pancy at long times is independent of the initial state of 
the phonon. Third, the term involving the summation 
over k in Eq. (j47|) decays to zero as well, resulting in a 
compact expression for the phononic occupancy at long 
times: 

ntit ^ ^) = A2 ^ l5(efc + {^k - ^of ■ (57) 

fc>0 

Finally, one can show that Eq. ([57)) is just the zero- 
temperature equilibrium phonon occupancy with respect 
to the full Hamiltonian,'^^ implying thermalization at 
long times. This result on its own is not surprising, 
since it has been rigorously shown by Ambegaokar^" that 
Hamiltonians involving a local bosonic mode coupled lin- 
early to a macroscopic bosonic bath do indeed equilibrate 
at long times in response to a local quantum quench. Be- 
low we analyze in detail the decay to the new thermal 
equilibrium. 

Figures [Tj and [2] summarize the time evolution of nb{t), 
for different coupling constants and different initial con- 
ditions. In Fig. [1] we have plotted nh(t) in response to an 
abrupt switching on of the electron-phonon coupling g, 
with the phonon initially occupying the empty state |0) 
at time t = 0. Different values of g are depicted. Starting 
from nf,(0) = 0, the time-dependent occupancy first over- 
shoots its new equilibrium value, to which it then decays 
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FIG. 1: (Color online) Time evolution of the phonon occu- 
pancy ni,{t) following an abrupt switching on of the electron- 
phonon coupling g at time t — 0, with the phonon initially 
occupying the empty state |0). Here ujo/Dd = 0.2, while g/T 
equals 0.229 (green), 0.28 (red), and 0.324 (black). The cor- 
responding values of (/■^/(oJoT) are 1/6, 1/4, and 1/3, respec- 
tively. Inset: A fit of the g/F = 0.324 curve to the functional 
form of Eq. (|58|) using the fitting range 9 < uJot < 95. The 
two curves practically coincide above ujot = 8. 

with superimposed oscillations. The osciUatory decay is 
well described by the long-time behaviors of Ii (t) , I2 (t) , 
and the term involving the sum over k in the expression 
for F{z,t). Indeed, based on our previous analysis one 
expects nb{t ^ ^/Dd) to follow the functional form 

nb{t) = [A sin(2m + 0) + S] e^^t/^o _^ /j^ (58) 

with Q, and tq equal to lu and r of Eqs. and ((55|) . The 
inset of Fig. ([T]) shows a typical fit of the g/T = 0.324 
curve to the functional form of Eq. ([55]) using the fitting 
range 9 < Lu^t < 95. While some deviations are seen 
at shorter times, the two curves are hardly distinguish- 
able above oJot = 8. Moreover, the extracted values of 
n/ujO = 0.896 and tqWo = 35.5 fall within 1.2% from 
those of oj and r quoted above. The agreement between 
the predicted and extracted parameters is equally good 
for the two curves with the smaller values of g, confirm- 
ing our analytic predictions for the long-time behavior of 
nb{t). 

Figure [2] displays the complementary dependence of 
nb{t) on the initial state of the localized phonon. As 
emphasized above, nb{t) depends on the initial state of 
the phonon via two parameters only: nb(0) and (66) (=0. 
Hence each curve with ribiO) > corresponds to a family 
of initial states. It is nevertheless useful to have a par- 
ticular initial state in mind by assigning a representative 
state to each combination of n;,(0) and {bb)t=o- One pos- 
sible choice of states for the two curves with nf,(0) > 
are [2|0) - i\2)]/V5 (red line) and |1) (green line). The 
curve with nb(0) = (black line) corresponds exclusively 
to |0) as the initial state. 

As in Fig.[Tl all curves in Fig.[2]can be fit equally well to 
the functional form of Eq. ([55)) using the same pair of val- 
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FIG. 2: (Color online) Time evolution of the phonon occu- 
pancy nb{t) following an abrupt switching on of the electron- 
phonon coupling g, for g/F = 0.324, oJo/Dd = 0.2, and differ- 
ent initial phononic states. Each of the curves with nb{0) > 
corresponds to a family of initial states whose values of nt{0) 
and {bb)t^o are specified in the legends. Representative states 
for each category are [2|0) - i\2)]/V5 (red) and |1) (green). 
The curve with n(,(0) = (black) corresponds exclusively to 
the initial state |0). 

ues for and tq that were extracted for ^^(0) = 0. Gen- 
erally speaking, the larger is nb{0) the more pronounced 
is the component of the pure exponential decay, while 
the magnitude of the superimposed oscillations is more 
sensitive to {bb)t=o- 

Another quantity of interest is the time evolution of the 
phonon displacement, Q{t). Since Q is strictly zero for 
Ed = in thermal equilibrium, its time evolution remains 
pinned to zero unless either or {b)t=o is nonzero. In 
this section we consider the latter possibility where {b)t=o 
is nonzero. A straightforward evaluation of Q(t) using 
Eq. ((33)) and its Hcrmitian conjugate yields 

Q{t) = V2Re{[h{t)+i;{t)]{b)t=o}, (59) 

whose dependence on the initial state is reduced to the 
sole parameter {b)t=o- Writing the latter in terms of 
its magnitude and phase, {b)t=o = \{b)\e^'^, the time- 
dependent displacement depends linearly on |(6) | . The 
dependence on ip is less transparent as it requires detailed 
knowledge of Ii{t) and l2{t)- Numerical calculations re- 
veal, however, that the dependence on ip is rather weak, 
hence we focus our attention hereafter on = 0. 

FigureOdepicts the time evolution of Q{t) for (6)t=o — 
1 and two representative values of the electron-phonon 
coupling g. As can be seen, Q{t) displays damped oscil- 
lations with a frequency and decay time that depend on 
the magnitude of g. Indeed, based on our previous anal- 
ysis of Ii (t) and I2 (t) one expects the long-time behavior 
of Q{t) to follow the functional form 

Q{t)=A sm{nt + (/i)e-*/^° (60) 

with f2 and tq equal to lo and r, respectively. Fits 
to Eq. (jSO)) yield excellent agreement, with values of 
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FIG. 3: (Color online) Time evolution of the phonon dis- 
placement Q{t), starting from an initial phonon state where 
{b)t=o = 1- Here uo/Dd equals 0.2. Two representative val- 
ues of the electron-phonon coupling are depicted: g/T = 0.324 
(black) and g/T = 0.229 (red), corresponding to (/^/(ljoT) — 
1/3 and 1/6, respectively. 



and To that coincide to within less than 1% with those 
extracted from Fig. [T] using fits to Eq. (|58)) ."^^ Thus, 
Q{t) displays a relaxation time twice as long as that 
of nb{t) and half the frequency of oscillations. Such a 
relation is quite natural for a classical oscillator where 
nb{t) = {b'^{t)){b{t)) - Q^{t), but is less obvious for the 
quantum case considered here. 



C. Phononic wave function 

Lastly we shall address the time evolution of the 
phononic wave function, defined as 



|V'ph(x,i)l' = ('5(Q(t)-.x)). 



(61) 



Here, Q{t) — W {t,0)QU{t,0) is the phonon displacement 
operator in its Heisenberg representation, a; is a dimen- 
sionless position coordinate, and averaging is taken with 
respect to the initial state of the system. In thermal 
equilibrium |'0ph(;E, was calculated by Doraj^^ who 
showed that it takes a simple Gaussian form. Below we 
extend the calculation to nonequilibrium quench dynam- 
ics, allowing for an arbitrary initial phonon state. 
Following Dora we begin by rewriting Eq. (|6ip as 



\ipph{x,t)\^ = 



(62) 



Since Q(t) is linear in bosonic operators, and since av- 
eraging on the right-hand side is taken with respect to 
a product state of the filled Fermi sea and the initial 
phonon state, (^e^^(Q(*)-^)) can be recast as the product 
of two independent averages of the conduction-electron 
and local-phonon components of Q{t). Explicitly, denot- 
ing the two components of Q{t) by Qc{t) and Qb{t) one 



JsQ{t)\ _ I isQb(t) 



)b(e 



FS : 



(63) 



where (. . .)fs a-nd (. . .)b stand for averaging with respect 
to the filled Fermi sea and the initial phononic state, 
respectively. 

Each of the two averages in Eq. (|63p can be evaluated 
in turn using standard bosonic techniques. Consider first 
the conduction-electron component. As the average is 
taken with respect to the ground state of a free bosonic 



bath one can use the identity {e^ 



,{A')12 



applicable 



to any operator A that is linear in bosonic creation and 
annihilation operators. This results in 



FS = exp 



FS 



(64) 



Moving on to the local phonon component, we first note 
that Qhii) has the explicit form 



(65) 



where 



h{t) - ^[/*(t)+/2(t)] =y2c.oA2 5]Cfcl5(efc + ^^)l' 
X [(cfe - wo)e"'=* + (efc + c^o)e-"'=*] . (66) 



Next we use the identity e'^"'"^ = e^e^ ey^^^l"^ ^ applica- 
ble to any two operators A and B whose commutator is 
a c-number, to write 

^g»sQi,(t)^j^ = g-(sV2)|/3(i)|'^gis73(t)''*gi'S-f3*(*)6)^. (67) 

The combination of Eqs. dHS]), dMl), and (HT]) then yields 



IV'ph(a^,i)r = 
with 



OO 7 

2 _ / ^-'y(t)s^-isx /^^isl3(t)b^ ^isn (t)b\ 

— OO 



27r 



b 

(68) 



7(0 = cj'oX' - + \\h{t)? (69) 



and 



K{z,t) 



k>Q 



q>0 



X 



-ie„t 



eq + Z 



(70) 



Equation allows one to calculate the phononic 
wave function for arbitrary time t > and any initial 
phononic state. Before turning to concrete examples let 
us address some generic features of |V'ph(a^, Since 
Islt) decays to zero as i — > oo, the expectation value on 
the right-hand side of Eq. (1551) reduces asymptotically to 
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FIG. 4: (Color online) Time evolution of the phononic wave 
function lipix, t)\'^, starting from an initial coherent state with 
A = 1. Here ujo/Dd = 0.2 and g/T ^ 0.324. 



FIG. 5: (Color online) Same as Fig. 21 starting from an initial 
phonon state where fib = 1. 



one regardless of the initial state of the phonon. Further- 
more, repeating the same type of analysis as beforehand 
one finds that the term involving the sum over q in the 
expression for K{ek — i^i^t) decays to zero with the re- 
laxation time r, and that "f{t) decays to its asymptotic 
value 
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lim J{t) = X^oJlJ2ek\9{ek + ^v)\' 



(71) 



k>0 



with the reduced relaxation time r/2. The phononic 
wave function thus takes the asymptotic Gaussian form 



|V'ph(a;)| = lim \tppi,{x,t)\^ 



1 



exp 



„2 1 



47 



(72) 



Lastly, recognizing that 7 is half the thermalized expec- 
tation value of Q^, i.e., 27 = ((5^)cq, we recover hereby 
the equilibrium result of Dora^^ at long times. 

While the asymptotic form of the phononic wave func- 
tion is independent of the initial state of the phonon, the 
associated relaxation time does depend on whether {b)t=o 
is zero or not. To see this we note that |'(/'ph(a;, i)P has 
two sources of time dependence originating from j{t) and 
Isit). While j{t) decays to its asymptotic value 7 with 
a relaxation time equal to t/2, l3{t) decays to zero with 
a relaxation time that is twice as long. The relaxation 
of |'0ph(a;,i)p depends then on whether the expectation 
value on the right-hand side of Eq. (j68| has a contribu- 
tion that is linear in l3{t) and I^it) or not. If {b)t=Q — 
there is no such linear contribution, hence |'0ph(a;,i)P 
approaches its asymptotic form with the relaxation time 
t/2. If, on the other hand, {b)t=o is nonzero then there 
is such a contribution and |?Aph(a;, i)p relaxes on a longer 
time scale equal to r. 

Although Eq. (1551) applies to any initial phonon state, 
of particular interest are those cases where the phonon 
initially occupies either a coherent state or an eigenstate 



of fih = b^b. If the initial state is a coherent state, i.e., 
b\^po) = X\ipo), then 



— e--^-'")b = e ^-3^";-/ ^ e-v^''^ (73) 



resulting in 



l^ph(a;,t)|' 



1 



exp 



47 (i) 



(74) 



The phononic wave function is therefore a simple Gaus- 
sian, characterized by the time-dependent average Q{t) 
and the time-dependent width a{t) = yj2^{t). If the ini- 
tial state is an eigenstate of fif, with the eigenvalue n the 
phononic wave function is somewhat more convoluted, 
given by the formal expression 



l^ph(x,i)l'= E 



2m j2m 



(75) 

Alternatively, Eq. (|75|) can be rewritten using Hermite 
polynomials as 



n 

mi m 



!V^[47(i)]™+i/2^"^2/je , 

(76) 

with y = x/^A-yit). 

The time evolution of the phononic wave function is 
displayed in Figs. U and [5] for two representative initial 
configurations of the phonon: a coherent state with A = 1 
(Fig. H} and an eigenstate of fn, with the eigenvalue n — 1 
(Fig. [5]). Starting from a coherent state, the phononic 
wave function evolves through a sequence of Gaussians, 
as can be seen in Fig. |4| The center of the Gaussian, 
Q{t), oscillates from at time i = to zero as t 00 
according to the black curve in Fig. |31 while its width 
increases form l/\/2 to 1. In contrast, the phononic wave 
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function undergoes a qualitative change in shape when 
starting from hj_ = 1. Here \'4'{x, t)p is initially composed 
of two symmetric peaks that gradually merge to a single 
Gaussian at long times. This behavior can be understood 
from the explicit form of Eq. (1761) with n — 1: 

(77) 

At t = one can show that \h{Q)\'^ = 27(0) = 1/2, 
resulting in 

|^ph(a;,t = 0)|2 = (78) 

As time increases |/3(t)p gradually decays to zero, leav- 
ing us with the thermalized Gaussian of Eq. (|72l) . The 
transition between the two forms of the wave function is 
therefore driven by the relaxation of |/3(i)p, which hap- 
pens on a time scale of r/2. 

D. Thermalization in the presence of integrability 

Tracking the time evolution of the phononic occupancy, 
displacement, and wave function in response to switching 
g on, we observed in the previous subsections that all 
quantities eventually approach their thermal equilibrium 
values with respect to the full Hamiltonian. In other 
words, the system thermalizes at long times. Indeed, 
this limit was rigorously shown by Ambegaokar^'' for a 
class of bosonic models that include our Hamiltonian of 
interest. More generally, it was shown by Doyon and 
Andrei^ in the context of interacting quantum dots that 

lim t/^(t,0)e-^^"C/(i,0) = e-^^, (79) 

f 00 

provided the bath is a large Fermi sea and % — "Ho is 
a local perturbation, as is the case here. However, one 
may wonder at this point how can an integrable system 
thermalize given the infinite set of conservation laws it 
possesses [specifically, the occupation numbers aj.afc, see 
Eq. (EHl)]? 

To address this question consider the quench dynamics 
starting from an excited Fermi sea 

/ J. \"fcl / 4. \ "few 

|»^/ci,---,rifc„)g=o = (^a[J 10)9=0, 

(80) 

obtained by creating several particle-hole excitations 
with momenta fci, ^2, . . . , fcAf above the filled Fermi sea of 
the unperturbed system. Starting from the product state 
IV'o) = Wb) (S" \nki-, . ■ . , nkr^)g=o of the excited Fermi sea 
and a local phonon state with the occupation number nt, 
one can repeat the calculation of Sec. IIVBI to track the 
time evolution of the phonon occupancy. At long time 
one finds 

nfc(<^oo) = X'^^£,'^\g{eq + ii])\'^ {eq-ujof 

q>0 



N 

= g(7iA:i,...,nfe„ |6'l'6|nfe,,...,nfe„)g, (81) 

where 

K.---,"fc«>9= (4J •••(«L) 10)9 (82) 

is the corresponding eigenstate of the full Hamiltonian, 
obtained by creating scattering-state excitations with 
identical quantum numbers above the ground state of 
the full system. Thus, while the initial state of the local 
phonon is wiped out in the course of the evolution, the 
quantum numbers characterizing the initial state of the 
Fermi sea are preserved. In other terms, the conserva- 
tion laws constrain the bulk but not the local degrees of 
freedom. Since scales as 1/L, local observables are 
independent of the initial state of the Fermi sea as long 
as the initial excitation energy is not extensive, i.e., 

1 ^ 

-Y^Uk^^O (83) 

in the thermodynamic limit. We therefore conclude that 
the local phonon thermalizes while the bath does not, 
and that the conservation laws which constrain the bath 
dynamics do allow for a generic evolution of local degrees 
of freedom. 



V. ABRUPT CHANGE OF PHONON 
FREQUENCY 

The second quench dynamics we consider is the re- 
sponse to a sudden change in the phonon frequency. 
Namely, the system is taken to occupy the ground state 
of the Hamiltonian of Eq. (|23|) at time t = when 
the phonon frequency is abruptly shifted from ujq to 
uji = LdQ + 6io > 0. In contrast to the electron-phonon 
coupling, which is difficult to control in actual devices, 
the frequency of vibrations can be tuned electrically in 
suspended carbon nanotubes4^ This offers a potential 
realization of the present scenario. For concreteness we 
restrict attention in this section to = and zero tem- 
perature, though both restrictions can be relaxed. Ac- 
cordingly, our interest will center on nh(t), as Q(t) is 
pinned by symmetry to zero. Similarly, the phononic 
wave function retains a Gaussian form centered about 
x = at arbitrary time t, with a time-dependent width 

equal to J{Q'^{t)). 



A. Time evolution of phononic operators 

The general strategy for calculating the time evolution 
of physical observables is similar to the one taken in the 
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previous section, except that the initial Hamiltonian is 
now the full Hamiltonian Ti of Eq. ([25]) with e^j set to 
zero, and the final Hamiltonian is given hy H' = H + 6% 
with 



(84) 



The technical details are slightly more cumbersome, 
though, since the time evolution of h^{t) is carried out 
by expanding in terms of the eigenmodes 7^ and 7^ of 
the final Hamiltonian H', whereas the evaluation of ex- 
pectation values requires an expansion of h^{t) in terms 
of the eigenmodes and of the initial Hamiltonian 
%. In other words, one needs to know how to convert 
from the eigenmodes of T-L' to those of T-L. 

There are two approaches one can take to achieve this 
goal. The first is to invert Eq. (PE)) and its Hermitian con- 
jugate in order to express a^, aj, b and b^ in terms of the 
eigenmodes of %, and to plug the resulting expressions 
into the expansion of 7]! in terms of a^, aj, b and b^ . An 

alternative approach is to directly express ^\ in terms of 
aq and by solving the modified Lippmann-Schwinger 
equation 



bin'] 



ir]{al 



(85) 



To this end, it is necessary to first write H' in terms of the 
eigenmodes of H, which follows directly from Eqs. ((29)) 
and (|32|) . As we prove in Appendix[3 the two methods of 
computation are equivalent, allowing us to use the latter 
approach which is more concise. Differing all details of 
the calculation to the Appendix we quote here only the 
end result: 
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al + 2SujX Ckgi^k + ill) 
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-g{eg + ir])- 



- WoWi 
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where 



5(^) 



2wiE(0) 



(86) 



(87) 



is the same function of Eq. (j27p with uq replaced by wi . 

Since H' has the same exact form as the Hamiltonian H 
of Eq. ([23]) only with ujq replaced by wi , one can borrow 
all results derived previously in Sec. IIIII for the latter 
Hamiltonian. In particular, Ji' is diagonal in the new 
basis set. 



k>0 



and the expansions of b^ and b'^{t) detailed in Eqs. ^ 
and p3p still hold upon substituting wi, g, 7^, and in 
for ojQ, g, ak, and a\, respectively. Plugging Eq. (1551) and 
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FIG. 6: (Color online) Time evolution of the phononic oc- 
cupancy following an abrupt shift in the phonon frequency 
from Ldo/Dd = 0.2 to lji = ljq + Sco with Slj = ±0.3a;Q. The 
electron-phonon interaction is held fixed at g/T = 0.324. 



its Hermitian conjugate into Eq. ([55]) we finally obtain 
the desired expansion of b^{t) in terms of the a^'s: 



b\t) = A^6[.g(efc-i77)(efc+cji)e*' 



k>0 



+2Sojgiek - iv)J{^k ~ t)] ct\ 



k>0 



+25ujg{ek + i'n)J{-ek - irj, t)] afe,(89) 
where we have defined the auxiliary function 

J{z,t) = A2^Cfel.9(efe + "7)l'[(efe+^i) 



k>0 



£k- z 



-h(efc-wi) ■ e 

ek + z 



B. Phononic occupancy 



(90) 



With Eq. ([89|) at hand, we are now in position to eval- 
uate expectation values pertaining to the local phonon 
mode 6^. Focusing on the time evolution of the phonon 
occupancy nb{t) = {b\t)b{t)), we note that {t)b{t) is 
quadratic in a\ and a^. Since the expectation value 
is taken with respect to the ground state of the ini- 
tial Hamiltonian T-L , the only nonzero contributions stem 
from the diagonal terms OLka\, resulting in 

ni{t) = A2^efcl.9(efc + i'7)(-efe+'^i)e-"** 

fc>0 

+25^g{ek + i^)J{~ek - «7?,t)|'. (91) 

At t = 0, Eq. (|9T|) properly reduces to the equilibrium 
expectation value of n;, with respect to % specified in 
Eq. ([57)) . To see this we note that Eq. ([M)) must coincide 
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at time t = with Eq. (I32l) . as both expressions offer 
an expansion of in terms of the a^'s and ct^'s. Equat- 
ing the corresponding expansion coefficients one finds the 
identity 

5(efe + ify)(-e/c + ^o) = '2Sujg{ek + iv)J(-fk ^ iv,^) 
+ 9{ek+iv){-^k+i^i), (92) 



from which the equivalence of Eq. (l57l) and Eq. ((9T|) 
at time t = immediately follows. In the opposite 
limit of long times Eq. (j9T|) reproduces the new equi- 
librium expectation value of hh with respect to H' . In- 
deed, using a similar analysis as beforehand one can show 
that J{—ek — if], t) decays to zero with the same relax- 
ation time T and frequency of oscillations w as listed in 
Eqs. ([51]) and ([55)1 . subject to the substitution of wq with 
LJi. This in turn leaves us at long times with 

^ oo) = ^ Cfc Vgi^k + i-nf (efc - ^i)', (93) 

fc>0 

which is the thermalized expectation value of fib with 
respect to T-L' . Thus, as expected, nb{i) interpolates be- 
tween the two equilibria expectation values of fib- 

Figure [S] shows the complete time evolution of nb{t) 
for two opposite shifts of the phonon frequency. Once 
again the curves take the form of damped oscillations 
with the relaxation time t/2 and the frequency of oscil- 
lations 2u). Consequently, the decay time and frequency 
of oscillations differ substantially between Sui — O.Swq 
and 5u = — O.Swq, in accordance with the substitution 
cjo ^ t^i = (1 ± 0.3)a;o in Eqs. ([Ml) and ([55]). The larger 
is (jJi the smaller are the new thermalized expectation 
value of hb and the amplitude of damped oscillations that 
nb{t) undergoes. 



VI. ABRUPT SHIFT OF ENERGY LEVEL 

The third and final quench scenario we consider is the 
response to an abrupt shift in the electronic energy level, 
which has been held fixed up until now at resonance with 
the Fermi energy. Specifically, we assume that the system 
resides at time t < in its ground state for e^j = 0, when a 
nonzero is suddenly switched on. This has the effect of 
breaking particle-hole symmetry, dynamically generating 
a nonzero displacement Q{t) of the local phonon along 
with deviations of the level occupancy from half filling 
[i.e., nd{t) 7^ 1/2]. These two observables will be our 
main focus of interest. Of the different quench scenarios 
under consideration the present one is by far the most 
accessible experimentally, as the energy level can be 
controlled quite efficiently using suitable gate voltages. 

The foundations for calculating Q{t) and nd{t) in this 
scenario have been laid down in Sec. IIIIBI Specifically, 
from Eq. (|38p one has that 



5(efe-*?7)(efe+wo)e"'^y 



k>Q 



-g{ek + ir])iek - wo)e "'"'/Sfc 



(94) 



which, when combined with Eqs. (|36l) and (l39l) . yields 



k>0 



-g(efc-f i?7)(efc-wo)e 



(95) 



with 



fc>0 



9 



(k + iri 
1 



Efc - if] 



(96) 



ttTwo 1 - 2g2/(7ra;or)' 

Accordingly, the phonon displacement is equal to Q{t) = 
\/2Re{/ii(t)}, which follows from the fact that at and aj. 
average to zero with respect to the initial state. Similarly 
from Eq. (PT|) one has that 



nd(t) 



1 



fc>0 



(97) 



which, when combined with Eqs. p6l) and (H^ . yields 



{t) = h2{t)+aJ2ik{el-u;^) \g{ek-iv)e'"'*al+R.c 



with 



h2{t) 



k>0 



(98) 



irT 1 - 2gy{T:LUoT) 



^daJ2^k\9i<^k+iv)\^i£l 



,2\2 



k>0 



f-k + iv ik-irj 



(99) 



Hence the occupancy of the localized level is simply given 
by Tidit) = h2{t). 

The occupancy nd{t) of the localized electronic level is 
depicted in Fig. [71 Several points are noteworthy. First, 
nd{t) — 1/2 depends linearly on ed in our solution. This 
property stems from the fact that couples linearly to 
the bosonic degrees of freedom, in accord with the as- 
sumption that \ed\ ^ r. Indeed, as \ed\ is increased the 
mapping onto the bosonic Hamiltonian of Eq. (j23p grad- 
ually breaks down, generating higher order corrections in 

Second, the dynamics of nd{t) is composed of two dis- 
tinct segments: fast dynamics on the scale of 1/Dd ~ 
1/r, where most of the charge redistribution takes place, 
followed by an extended region of damped oscillations. 
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FIG. 7: (Color online) Time evolution of the occupancy nd{t) 
of the localized electronic level, following an abrupt change 
in its energy from ed = to ^ 0. Here uo/Dd = 0.2 and 
g/r — 0.324. Note that nd{t) — 1/2 depends linearly on e^- 
The green dashed line shows a fit to Eq. (llOip using the fitting 
range 10 < ujot < 100. Inset: A zoom in on the short-time 
behavior. The red dashed curve shows the analytical form of 
Eq. (fTOO|l . 



The short-time dynamics originates from the high-energy 
end of the summation over k in Eq. (1991) , and is given by 



^ "2 Trr 1 + {tDdy 



(100) 



(see inset of Fig. [T]). 
from the exponential high-energy cutoff imposed by ^j. 



This simple analytical form stems 

2 
k 

[see Eq. (12^ ]. While the functional form of nd{t) may 
differ from Eq. (|100p for other cutoff schemes, the rele- 
vant time scale t ~ l/F and the characteristic change in 
occupancy Srid e^/Trr experienced within this time seg- 
ment should be generic. As for the damped oscillations, 
these are expected to take the functional form 



nd{t) — A sm{Qt 



C 



(101) 



with fl and tq equal to uj and r. We confirm this form 
in Fig. [3 where the fitted values of and tq agree to 
within 0.1% with those extracted from Fig. [T]by fitting 
nb{t) to Eq. dSHl.'^i 

Lastly, from the asymptotic long-time behavior of nb{t) 
one can deduce the dimensionless parameter controlling 
the perturbative expansion in g in thermal equilibrium. 
For g = and arbitrary e^, the exact equilibrium oc- 
cupancy of the electronic level is given by the standard 
expression 



1 1 /Cd 
nrf = arctan — 

2 TT VF 



(102) 



which reduces to = 1/2 — e^/TrF for \ed\ ^ F. This 
latter result is accurately reproduced by our treatment 
upon setting g equal to zero. For nonzero g we find that 



1 



7rri-252/(^^or)' 



(103) 



FIG. 8: (Color online) Time evolution of the displacement 
Q{t), following an abrupt change in the level energy from 
Ed = to 7^ 0. All model parameters are the same as in 
Fig. El For comparison, the red curve plots Q{t) of Eq. p06|l . 



revealing that the true expansion parameter is g'^ /luoT, 
i.e., the ratio of the polaronic shift g^ /uiq to the hybridiza- 
tion width F. The effect of the electron-phonon coupling 
g is to increase the deviation from half filling for a given 
value of edi signaling a narrowing of the electronic reso- 
nance according to 



F -> Feff = F 



TTU/Q ' 



(104) 



The time evolution of the phonon displacement is plot- 
ted in turn in Fig. [5] Similar to the level occupancy, Q{t) 
depends linearly on and undergoes damped oscillations 
with the relaxation time r and frequency lo. It lacks, 
however, the fast dynamics that the level occupancy ex- 
periences on the time scale of 1/F. In equilibrium rid and 
Q are related through 



1 

Ud-- 



(105) 



which is an exact result applicable to arbitrary e^, wq, 
and g.'^* It is thus natural to ask whether this general 
relation extends to nonequilibrium dynamics. To this 
end, in Fig. |S] we have plotted 



Q(t) = -V2 



9_ 

OJQ 



nd{t) 



(106) 



alongside Q(t). Although nearly in phase, the two quan- 
tities are characterized by vastly different amplitudes of 
oscillations, marking the breakdown of Eq. (jl05p under 
nonequilibrium dynamics. The latter relation is restored 
only asymptotically as the system thermalizes. 



VII. DRIVEN DYNAMICS 

Up until now we considered the response to a single 
quantum quench. Quantum control of nanodevices often 
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requires the usage of driven dynamics, where periodic 
forcing is applied to the system. Such drives are a the- 
oretical challenge to describe since the system not only 
remains permanently remote from thermal equilibrium, 
but it never even reaches steady state. Remarkably, we 
are able to extend our exact solution to a rather broad 
class of driven dynamics where the forcing couples lin- 
early to the bosonic degrees of freedom. As we discuss 
below, this class of drives includes at least two physically 
relevant scenarios where periodic forcing is applied either 
to the localized phonon or to the electronic level. Accord- 
ingly, we begin our derivation with a general discussion 
of this class of drives before turning to the two concrete 
examples of interest. 



A. Drives that couple linearly to bosons 

The general setting we consider consists of a system 
that resides at time t < in thermal equilibrium, when a 
time-dependent drive is suddenly applied to it. In formal 
terms, the Hamiltonian of the system is changed abruptly 
at time < = from the Hamiltonian H of Eq. (|23p to 
n'it) =-H + -Hdrivc(t) with 



^drivc (t) = ^[Mfe(i)al-fM,!(tK 

fc>0 



(107) 



Here we have assumed that the drive couples linearly to 
the eigenmodes of H, exploiting the fact that any lin- 
ear combination of the original bosonic degrees of free- 
dom can be expanded in terms of the scattering-state 
operators ak and a\. The coefficients Mk{t) depend on 
the exact scenario under consideration and will typically 
have the separable form Mk{t) — A{t)mk- Nevertheless, 
we shall regard them for the time being as general co- 
efficients without making any further assumption about 
their form. The initial value of the energy level will 
be taken for simplicity to be zero, though the extension 
to nonzero is quite straightforward. 

A convenient way to incorporate the time-dependent 
drive is via the Heisenberg equation of motion for the 
scattering-state operators, which takes the form 



ctkit) = -iekUkit) - iMk(t), 



(108) 



subject to the initial condition ak{t = 0) = ak- Here the 
ffist term on the right-hand side of Eq. (jl08l) is due to 
T-L and the second term is due to 'Hdrivo- Equation (jlOSp 
has the formal solution 



lie 

10 



'"'^''-''>Mk{t')dt', (109) 



from which the time evolution of all physical operators 
of interest can be deduced. For example, combining 
Eq. ([5^ with Eq. (|109p and its Hermitian conjugate one 
obtains 



b^{t)^bl{t) + iXBit), 



where fcj {t) is the time-evolved operator in the absence of 
a drive [see Eq. p3p ] and B{t) is a time-dependent shift 
given by 



Bit) = J2^k 



k>0 



g{ek^iv){ek+ooo) / A4!(t')e-'"'=(*'-*)di' 



,(111) 



+g{ek + iv){ek~ujo) / Mk{t')e'"'^''~''>dt' 



Accordingly, the phonon displacement takes the form 

Q{t) = -V2\lm{B{t)}, (112) 



while its occupancy reads 



(113) 



Here ni^^ denotes the equilibrium phononic occupancy in 



the absence of a drive, given by Eq. ([57| for T = 0. Note 
that in deriving Eqs. (|112p and ()113p we have made use 
of the fact that b{t) + b'<{t) and b''{t)b{t) are averaged with 
respect to the equilibrium density operator corresponding 
to H, which is diagonal in the occupation numbers a\ak- 
As a result (ak) and (aj.) identically vanish. A similar 
calculation for the time-dependent occupancy dnd(t) = 
nd{t) — 1/2 of the localized level yields 



5nd = a^Cfc(efc-t^o)[-*5(efc + i?7) / Mfc(i') 
fc>o 



+ ig{ek-iv) / Ml{i!)e'^'-^''-'Ut']. (114) 
Jo 



Equations (I113p and (|114p combined provide us with 
a formal solution for the occupancies of the local phonon 
and the electronic level for a general drive. Furthermore, 
these expressions apply to arbitrary temperature T, pro- 



ided 



,(0) 



is taken to be the equilibrium phononic oc- 



cupancy at that temperature. Below we utilize these ex- 
pressions to analyze two cases of practical interest, where 
ac forcing is applied either to the local phonon or to the 
localized level. 



B. ac forcing of the local phonon 

In the first scenario to be analyzed, ac forcing is applied 
at time t > Q to the local phonon, as described by the 
Hamiltonian term 



,{t) ^ Asin(m)(&^ +5). 



(115) 



(110) 



Here A, which has dimensions of energy, is the amplitude 
of the drive and Vt is the driving frequency (not to be 
confused with the fitting parameter previously used for 
analyzing the damped oscillations). Such forcing can be 
applied, e.g., to polar molecules using an ac electric field. 

Using the expansion of &^ in terms of the scattering- 
state operators given in Eq. p2p . the Hamiltonian term 
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of Eq. (IllSp can be recast in the form of Eq. (|107|) with 
the coefficients 



Mkit) ^ 2woAAsin(m)Cfcg(efe - ir]), 



(116) 



such that 



M*(i')e-'^^(*'-*'dt' = woAAaff(efc + tvKkit) (117) 



with 



a(i) 



efe - 17 



efe + 



(118) 



For computational convenience it is useful to add an 
infinitesimal imaginary part irj to the denominators in 
Eq. (|118p . thereby rewriting (k{t) as^ 



Ck{t) 



efc — + irj efe + + irj 
A somewhat lengthy calculation then gives 
A 



(119) 



[F{--n-ir],t)-F{n-ir],t)], (120) 



XB{t) 



where F{z,t) is the same function defined in Eq. (l47t . 
The phonon displacement and occupancy therefore take 
the rather compact forms 



Q{t) 



V2 



lu\{F{-Q.-iri,t)~ F{n~i7^,t)} (121) 



and 



A2 

— \F{-Vt-i7^,t)~F{Vt- 



-Z77,i)p. (122) 



The general structure of Q{t) and nb{t) could be un- 
derstood from properties of the function F{z,t). Since 
F{z, 0) identically vanishes for arbitrary z, the phononic 
occupancy and displacement properly reduce at time 
i = to their thermal equilibrium values, as they phys- 
ically should. As soon as t > 0, the two components 
of F{±Vt — irj) behave markedly differently. The first 
term in Eq. (j47p oscillates indefinitely with frequency 
n, whereas the term involving the sum over k under- 
goes damped oscillations with the relaxation time r and 
frequency ui of Eqs. (|55)) and (|54)) . respectively. Thus, 
there is a clear distinction between the roles of the two 
terms: while the first term in Eq. (j47l) survives at long 
times and is responsible for the long-time behavior, the 
second term contains all transients that decay in time. 
This leads to the following characterization of Q{t) and 
nb(t). At short times, t < t, the phonon displacement 
comprises of two components oscillating at frequencies 
and oj, while nb{t) contains four distinct oscillatory terms 
with the frequencies 2f2, f2 ± w, and 2uj. At long times, 
T ^ t, the phonon displacement oscillates with frequency 
f2 about zero, while ni,{t) oscillates with frequency 2f2 



C^ 60 
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FIG. 9: (Color online) Time evolution of 5nb{t) = ni,{t) — n^"' 
in response to ac forcing of the phonon according to 'Hdrivc(i) 
of Eq. (fTTSj) . Here ujo/Da = 0.2, g/T = 0.28, and fl = l.Sujo. 
Initially there is a rich structure involving the interference of 
four distinct frequencies. As transients decay (on a time scale 
of r), Sni,{t) gradually reduces to a single harmonic with a 
frequency of oscillations equal to 2Q. Inset: Zoom in on the 
earlier time segment uJot < 50, including a comparison to the 
stronger coupling strength g/T = 0.324 (red curve). 



about a new time- averaged value n;,. Explicitly, Q{t) 



and Sni,(t) = nf,(i) — reduce at long times to 



Q{t) ^ V2AiOQ\g{n + iT])\sm{nt - 



(123) 



and 



Sribit) ^—\g{n+ir])\^ [ujl+fl^ + {fl^ - cj^) cos(2m+20)] , 

(124) 

with (f) = arg{f/(r2 — ir])}. 

On physical grounds one expects the response to an ac 
drive to reduce at long times to a periodic function of 
time, containing all harmonics of the driving frequency 
fl. Surprisingly, the oscillatory parts of Q{t) and Snbit) 
consist in this limit of just a single harmonic each. While 
Q{t) tracks the driving field with a phase difference of (j), 
Snb{t) oscillates with the doubled frequency 2f2, lacking 
any signal at the principal harmonic fl. Such behavior 
is quite atypical, as is the absence of higher harmonics. 
We expect both features to qualitatively change as the 
electron-phonon coupling is increased beyond the validity 
of our solution. We further note that the doubling of 
frequency in Snb{t) is reminiscent of a similar doubling of 
frequency in the damped oscillations that ribit) undergoes 
in response to a quantum quench (see, e.g.. Figs. [T] and 
[3] and their accompanying texts) . 

Focusing on Snhit), its amplitude depends in a simple 
quadratic manner on A. Other than setting the overall 
amplitude, A has no additional effect on Snb{t). By con- 
trast, the shape of Snb{t) is quite sensitive to the driving 
frequency il, as demonstrated in Figs. IH] and [TUl When 
fl is tuned off-resonance with the intrinsic frequency uj 
of the system, see Fig. [SI the transient behavior shows a 
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FIG. 10: (Color online) Same as Fig. (9] with Q tuned to 
the resonance frequency: uu — 0.896a;o for g/T = 0.324 and 
(jj = 0.923a;o for g/V = 0.28. All other model parameters are 
the same as in Fig. [9] Note the vastly different vertical scale 
as compared to that used in Fig. O 

rather rich structure that stems from the interference of 
the four underlying frequencies 2D,, D,±uj, and 2uj. Only 
after all transients have decayed on a time scale of r does 
Snf,{t) approach its asymptotic long-time form of a single 
harmonic with the doubled frequency 2fl. 

A rather different picture is recovered when fl is tuned 
to the resonance frequency w, see Fig. [TOl Here both 
the transient and long-time behaviors are governed by 
the same single frequency 2fl = 2a;, resulting in much 
smoother curves. Quite striking is the substantial in- 
crease of the amplitude of oscillations upon approaching 
the resonance frequency. Indeed, the amplitude of the 
long-time oscillations is roughly two orders of magnitude 
larger in Fig. [TU] as compared to Fig. [HJ which is read- 
ily understood from Eq. (|124p . Since the amplitude of 
oscillations is given at long times by 

At^^\g{n + zrj)\'\u;l-n% (125) 

it displays a sharp resonance for « w where |(7(r2 -|- 
irj) I is sharply peaked.^^ The amplitude of oscillations at 
resonance can be crudely estimated as 

^r~^l-o'--^l, (126) 

with to and r approximately given by Eqs. (|54p and (1551) . 
respectively. Another interesting observation is the van- 
ishing of Ab for 17 = Wo. A plot of the amplitude Ab of 
the long-time oscillations as a function of is depicted 
in Fig. [m 



C. ac forcing of local electronic level 

In the second scenario that we analyze, ac forcing is 
applied at time i > to the localized electronic level, as 



FIG. 11: (Color online) The asymptotic long-time amplitude 
of oscillations At vs the driving frequency Q, for uJo/Dd = 0.2 
and two different strengths of the electron-phonon coupling 
g. Inset: A zoom in on the vicinity of the resonance peak. A 
logarithmic scale is used for the y ordinate so as to emphasize 
the vanishing of Ab for Q. = ujq. 

described by the Hamiltonian term 

■Hdrivc(t) = Asin(m) [rid - ^ ■ (127) 

Here, as before, A denotes the amplitude of the drive and 
is the forcing frequency. Experimentally such a drive 
can be realized by applying microwave voltage to a near- 
by plunger gate, similar to the setups used by Elzerman 
et al.'^- and by Kogan et al.— in their respective studies 
of the ac Kondo effect in semiconductor quantum dots. 

Using the mode expansion of Eq. ([M)) one can again 
recast the Hamiltonian term of Eq. (|127p in the form of 
Eq. (|107p . this time with the coefficients 

Mk{t) = Aasm{nt)^k9{ek ~ iv){4 ~~ ^o)- (128) 

Repeating the same sequence of steps detailed in 
Eqs. (|117p - (|119p and plugging the resulting expressions 
into Eq. (|114p . one obtains after a rather lengthy calcu- 
lation 

Snait) = nait) - i = ^lm{F{n - ir], t)} (129) 
2 5 

with 

F{z,t) = {z'-UJl)g{zMz)e'^' + X'Y.^'\9{ek + ^v)\' 

k>0 

x(z2-Wo2)M^ + ^ . (130) 

The function F(z,t) has similar properties to those of 
F(z,t). At i = it vanishes identically for any value 
of z, and is composed of two distinct components for 
t > 0: one that oscillates indefinitely with frequency 
ri, and another that oscillates with frequency w and de- 
cays with the relaxation time t for any z in the lower 
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i2/(n„ = 0.896 




FIG. 12: (Color online) Time evolution of Snd{t) in response 
to ac forcing of the local electronic level according to 'Hdrivc(t) 
of Eq. ([mil . Here uJo/Dd = 0.2 and g/T = 0.324. Two 
driving frequencies are shown, one {Q. = I.Scjq) off resonance 
and the other (Q. — 0.896aJo) on resonance with the internal 
frequency lu. 

half plane. Since 6nd{t) is proportional to the imaginary 
part of F{n — irjjt) [as opposed to 5ni,{t) that depends 
quadratically on F{±D,—ir], t)] it comprises at short times 
t < r of two oscillatory terms, one with frequency fl and 
another with frequency uj. As t exceeds r the latter com- 
ponent is progressively suppressed and Snd{t) gradually 
approaches its asymptotic long-time form 

Snd{t) = Adsm{nt + ip), (131) 

with 

Ad^^\{Q^- uDgin + i7])i:{Q + Z77)| (132) 
9 

and (fi = arg {(fi'^ — ujQ)g{il — iri)Y.{il — irj)^ . As before- 
hand, the long-time oscillations develop a resonance for 
fa UJ, albeit with a reduced amplitude as compared to 
Ab of Eq. (|125p . This reduction in amplitude stems from 
a weaker linear dependence of Ad on \g{Q + ir])\. Similar 
to Ab, the amplitude of oscillations is suppressed to zero 
for Q = ujo, leaving no signal at long times for this par- 
ticular frequency. A summary of our results is presented 
in Figs. [H and [H 

VIII. SUMMARY AND CONCLUSIONS 

In this paper we have presented an asymptotically 
exact solution for the nonequilibrium dynamics of a 
single-molecule transistor in response to various quan- 
tum quenches and drives. Our solution, which is based 
on a controlled mapping of the original Hamiltonian of 
Eq. ^ onto a form quadratic in bosonic operators,^'' 
is formally confined to weak electron-phonon coupling 
and near-resonance conditions for the electronic level: 
r >• ma.x{g,\ed\,g'^/uJo}. While some aspects of this 




FIG. 13: (Color online) The asymptotic long-time amplitude 
of oscillations Ad vs the driving frequency SI, for ujo / Dd = 0.2 
and two different strengths of the electron-phonon coupling g. 



regime can be accessed using ordinary perturbation the- 
ory in g, the ability to sum all orders exactly allowed 
us to (i) explicitly show how the system thermalizes fol- 
lowing a quantum quench, (ii) identify the different time 
scales that govern the dynamics of the system, and (iii) 
access the asymptotic long-time response to a periodic 
drive. 

Transient behaviors following a quantum quench were 
found to involve two characteristic scales^i — an intrinsic 
frequency uj and a relaxation time r approximately given 
by Eqs. ([M)) and (f55|) . respectively. Quite surprisingly, 
some observables, such as the phonon displacement and 
the electronic occupancy of the localized level, display 
damped oscillations with frequency w and the relaxation 
time T, while other observables, such as the phononic oc- 
cupancy, oscillate with frequency 2uj and decay with the 
reduced relaxation time t/2. The distinction has to do 
with the bosonic representation of the observable in ques- 
tion. If the latter is expressed as a linear combination of 
bosonic operators, the relevant frequency and decay time 
are lo and r. If, on the other hand, the the observable in 
question is quadratic in bosonic operators, the relevant 
frequency and decay time are 2a; and t/2, respectively. 

A special feature of our solution is the nature of the 
long-time response of observables to ac drives, whose os- 
cillatory component reduces to just a single harmonic 
with an amplitude that depends in a simple power-law 
fashion on the forcing amplitude A. The absence of 
additional harmonics in the long-time ac response is a 
direct consequence of the mapping onto a free bosonic 
Hamiltonian with a forcing field that couples linearly to 
the bosonic modes. Physically this implies that other 
harmonics, which are generally expected to exist for the 
original Hamiltonian of Eq. ([T]), are parametrically small 
for r ^ ra.a3i{g , g^ / . As the electron-phonon inter- 
action is increased such that max{g,g^/wo} approaches 
r, additional harmonics are expected to gain significance, 
that is provided the forcing amplitude A is not too small. 
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Concomitantly, the amplitudes of the different harmonics 
should gradually acquire a more elaborate dependence on 
A beyond a simple power-law form. We emphasize that 
this regime can no longer be described by the bosonic 
Hamiltonian of Eq. ((23)) . 

It would be interesting to compare our results with 
a numerical evaluation of the quench dynamics us- 
ing, e.g., the time-dependent numerical renormalization 
group (TD-NRG)i^ii^ Since lo and r are independent of 
the high-energy cutoff used in the electronic Hamiltonian 
of Eqs. ([TJ and (0), this should facilitate a direct compar- 
ison between the two approaches on time scales exceeding 
1/Dd ~ 1/r. The precise forms of w and r, as well as 
the short-time dynamics up to t ~ l/F, do depend on 
the cutoff scheme used for the bosonized Hamiltonian of 
Eq. Nevertheless, we expect our weak-coupling ex- 

pressions for w and t to apply in their present forms, 
as these coincide with low-order perturbation theory in 
g when applied directly to the electronic Hamiltonian of 
Eq. ((B). 

The true power of the TD-NRG lies, however, in 
its ability to treat arbitrary couplings strengths, which 
should enable one to go beyond the weak-coupling regime 
covered in this paper. It would be particularly interest- 
ing to see which aspects of our solution persist away from 
weak coupling, and what are the new qualitative features 
that are introduced as the electron-phonon coupling is in- 
creased. The study of stronger couplings along these lines 
is left for future work. 



1. A level at resonance with the Fermi energy 

We begin with a level at resonance with the Fermi 
energy, corresponding to the Hamiltonian of Eq. 
with td = 0. Our objective is to solve the Lippmann- 
Schwinger equation 



(Al) 



where 77 — >■ 0^ is a positive infinitesimal. To this end, 
we employ the methodology developed in Ref. Issl . Intro- 
ducing the Liouville operator CO = [0,7^], Eq. (|A1[) is 
rewritten in the form 



(£ + Efc + irf)ot\ ^ 
which has the formal solution 



(A2) 



irj 



C + ek + irj " 
Next we divide the Hamiltonian H into three parts 



(A3) 



k>0 

q>0 



(A4) 

(A5) 
(A6) 



and associate each Hamiltonian term with its own Liou- 
ville operator: £„0 = [0,'H„] (n = 0,1,2). Using the 
operator identity 
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C + ek+iri 
in combination with 

(A + ek 



+ fife + 47? 



- {C, + £2 



Cial = 0, 



1 



Co 



ek + ti] 
(A7) 



(A8) 



(A9) 



and 



Appendix A: Solution of the scattering-state 
operators 



In this Appendix, we detail the solution of the 
scattering-state operators for the different cases covered 
in the main text. Altogether three cases are consid- 
ered: (i) a level at resonance with the Fermi energy, i.e., 
ed = 0, (ii) a level off-resonance with the Fermi level, 
i.e., ed 7^ 0, and (iii) a local phonon with the shifted 
frequency uji ^ luq + 6uj. As emphasized in the main 
text, in the latter case we are interested in expanding 
the scattering-state operators corresponding to the fre- 
quency LOi in terms of those corresponding to the original 
frequency loq. 



C2al = -X^k{b^ +b), 
Eq. (jA3p is recast in the form 

1 



(6U6). 



(AlO) 



' C + ek+ir] 

Equation (jAllI) features two unknown quantities. 



(All) 



A, 



1 



C + ek + ir] 



-b^ and Bi. 



1 



C + ek+irj 



b. (A12) 



Our next goal is to explicitly compute these two operators 
by expressing them as the solution of two coupled linear 
equations. Once at hand, the scattering-state operator is 
simply given by al ^ al + X£,k{Ak + Bk). 
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To find Ak and Bk we resort once again to the op- 
erator identity of Eq. (IA7|) . Carrying out the relevant 
commutators one obtains the pair of equations 

{ek - i^Q + iv)Ak = + A VCgT— ^— -— (flg + ag), 

t^o ^ + 'k + iV 

(A13) 

{ck +UJO + iv)Bk = b - A ^ ^ {al + Uq). 

Z-/ ~r I 



q>0 



(A14) 

Applying yet again the operator identity of Eq. (IA7p to 
the right-most term in Eqs. (jA13[) and (jA14[) one arrives 
at 

{ek-iOa + iv)Ak = + ^{ek + iv){Ak + Bk) + Ck, (A15) 



{ek+ujo + iv)Bk = b-j:{ek + iv){Ak + Bk)-Ck, (A16) 
where 

^„ I 2; - e« z + e„ 

q>Q 



(A17) 



is the phononic self-energy and Ck equals 



— A 



9>0 



Cfc - eq + ii] ek + £q + if] 



(A18) 



Here in deriving Eqs. (jAlSp and (jA16l) we made use of 
the fact that 



1 



1 = 



1 



it 



(A19) 



Co + ek + iri \ o-q J Cfe T e<j + V "9 
Finally, introducing the 2x2 phononic Green function 

— S(z) — z — ojQ — S(z) 
Eqs. (lAlSp and (jA16p are rewritten in the compact form 



G(z) 



(A20) 



b^ + Ck 
b~Ck 



whose solution is 



(A21) 



(A22) 



Here CTz is the Pauli matrix. The scattering-state op- 
erators specified in Eq. ([26]) are obtained by combining 
Eqs. (UTTj) . (|H2)) . and (|X22)) . Note that the function 
g{z) defined in Eq. (j27p is simply minus the determinant 
of G-i(2). 



2. Extension to nonzero ed 

The case of a level ofi-resonance with the Fermi energy 
can, in principle, be treated using the same machinery as 
the one employed for = 0. We, however, shall present a 
more concise derivation that makes use of the scattering- 
state operators obtained for e^; = 0. As in the main text, 
the notation a\ will be reserved for the scattering-state 
operators when = while the new operators for 7^ 
are denoted by 

Our starting point is the formal solution 



Pi 



C + ek+irj 



(A23) 



where C pertains this time to the full Hamiltonian of 
Eq. ([231) with id ^ 0. Using the notations of Eqs. (R4)) - 
(|A6p . we divide the full Hamiltonian into two parts: 
■He .=0 ='Hq + 'Hi +'H2 and 



g>0 



q r 



(A24) 



Denoting the corresponding Liouville operators by Ce^=o 
and Ce^ , respectively, we employ the operator identity 



1 



1 



1 



C + Sk + 11] 



C + Sk+ir] 
to rewrite Eq. (IA24|) in the form 
1 



1 



£ed=o + ek+ir] 
(A25) 



Pi 



1 



-A 



C + ek+iv~'''' 
Recognizing that 

IT] 



IT] 



Cea=a + ek + iv 



-al (A26) 



we thus arrive at 



4 = 4' 



1 



C + €k+iri 



(A27) 



(A28) 



Since a\. and are both linear in the original bosonic 
degrees of freedom, their commutator is a simple c- 
number: 



[1 + gifk + ir])2ujQT,{ek + if])] 
-id^kg{ek + iv) (4 - ^0) • (A29) 



Consequently, using Eq. (jA28p . 

ek + it] 
which is precisely Eq. (1361) . 



(A30) 
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3. Change in frequency from to iui = iuo + Suj 

Lastly, we wish to expand the scattering-state opera- 
tors 7^ corresponding to the Hamiltonian H' = H + SH 
in terms of those corresponding to H alone (i.e., the at's 
and a\'s derived in Sec. lA Here H is the full Hamil- 
tonian of Eq. (1231) with id ~ and 6H equals 



Sn = SLub%. 



(A31) 



As in the previous subsection, we begin from the formal 
solution 



IT] 



+ ek +iri 



4, 



where C is the Liouville operator associated with H' . 
Denoting the Liouville operators corresponding to H and 
SH by £ and £slu , respectively, we make use of the oper- 
ator identity 



1 



1 - 



1 



+ Cfc + ill 



C + ek + iv 
to rewrite Eq. (jA32|) in the form 
1 



1 



C + ek 



1 



£ 



£' + ek+ irj 
Recognizing once again that 
ir/ 



17] 



£ 



+ Cfe + if] 



£ + ek + if] 



we arrive at 



7i 



1 



-£ 



£' + ek + IT] 



Suj^k ' 



- 11] 
(A33) 



(A34) 



(A35) 



(A36) 



which is analogous to Eq. (|A28[) of the previous subsec- 
tion. 

It is straightforward to confirm that Eq. (IA36[) is equiv- 
alent to the modified Lippmann-Schwinger equation 



bin'] 



(A37) 



Its usefulness stems from the fact that it allows one to 
directly expand ^\ in terms of the a^'s and aj's without 



resorting to the separate expansions of ^\ and a\ in terms 
of the original bosonic degrees of freedom. Indeed, using 
Eq. ([32]) and its Hermitian conjugate one has that 

£5loOi\ = 6uj\£_k9{<^k + irf) [(efc - '^o)^ - (efc + wo)6^] , 

(A38) 

such that 



with 



(A32) and 



A'k 



Bi. 



1 



£' + Ek + ill 
1 

£' + Efc + 



(6t - b) 
-I- h) . 



(A39) 



(A40) 



(A41) 



Similar to the derivation in Sec. lA 11 A'^ and B'^. are 
computed by expressing them as the solution of two cou- 
pled linear equations, obtained by applying the operator 
identity of Eq. (|X33)) to each of Eqs. (1X301) and (jUTI) . 
After some lengthy but straightforward algebra one ob- 
tains 



M{ek + iri) 



A' 
^k 

Bi 



fJ-k 



with 



M{z) = 



LOq 



{z^g{z) - 1) -Sujzg{z) 
—Sujzg(z) 1 — Sujujogiz) 



(A42) 



(A43) 



^fe 2A ^ £,geq 

q>0 



ek~ eq + irj ' ffe + eg + ir] 



and 



(A44) 



Vk = 2woA ^ 

9>0 



g{eq-ir]) ^ g{eg + iij) 

—a' H —an 

ek - eq + ir] ^ ek + eg + ir] 



(A45) 

Inverting the matrix M(ek + if]) to extract A'^, and 
and substituting the resulting expressions into Eq. (jA39l) , 
one recovers Eq. ([86| for 7^. 
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